Full configuration interaction approach to the few-electron problem in artificial atoms 
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We present a new high-performance configuration interaction code optimally designed for the 
calculation of the lowest energy eigenstates of a few electrons in semiconductor quantum dots (also 
called artificial atoms) in the strong interaction regime. The implementation relies on a single- 
particle representation, but it is independent of the choice of the single-particle basis and, therefore, 
of the details of the device and configuration of external fields. Assuming no truncation of the Fock 
space of Slater determinants generated from the chosen single-particle basis, the code may tackle 
regimes where Coulomb interaction very effectively mixes many determinants. Typical strongly 
correlated systems lead to very large diagonalization problems; in our implementation, the secular 
equation is reduced to its minimal rank by exploiting the symmetry of the effective-mass interacting 
Hamiltonian, including square total spin. The resulting Hamiltonian is diagonalized via parallel 
implementation of the Lanczos algorithm. The code gives access to both wave functions and energies 
of first excited states. Excellent code scalability in a parallel environment is demonstrated; accuracy 
is tested for the case of up to eight electrons confined in a two-dimensional harmonic trap as the 
density is progressively diluted up to the Wigner regime, where correlations become dominant. 
Comparison with previous Quantum Monte Carlo simulations in the Wigner regime demonstrates 
power and fiexibility of the method. 

Keywords: configuration interaction; correlated electrons; artificial atoms; Wigner crystal 



I. INTRODUCTION 

Semiconductor quantum dotsiiSi^ (QDs) arc nanome- 
ter sized regions of space where free carriers are confined 
by electrostatic fields. Typically, the confining field may 
be provided by compositional design (e.g., the QD is 
formed by a small gap material embedded in a larger- 
gap matrix) or by gating an underlying two-dimensional 
(2D) electron gas. Different techniques lead to nanome- 
ter size QDs with different shapes and strengths of the 
confinement. As the typical de Broglie wavelength in 
semiconductors is of the order of 10 nm, nanometer con- 
finement leads to a discrete energy spectrum, with energy 
splittings ranging from fractions to several tens of meV. 

The similarity between semiconductor QDs and nat- 
ural atoms, ensuing from the discreteness of the 
energy spectrum, is often pointed outi^iSiLSiS Shell 
structurQ£i2ii£ and correlation effectaSiiiiiS are among 
the most striking experimental demonstrations. More 
recently, the interest in new classes of devices has fu- 
elled investigations of QDs which are coupled by quantum 
tunneling^ In addition to their potential for promising 
applications, these systems extend the analogy between 
natural and artificial atoms to the molecular realm (for 
a review of both theoretical and experimental work see 
Ref.lll. 

Maybe the most prominent feature of artificial atoms is 
the possibility of a fine control of a variety of parameters 
in the laboratory. The nature of ground and excited few- 
electron states has been shown to vary with artificially 
tunable quantities such as confinement potential, den- 
sity, magnetic field, inter-dot coupHngi ^i^i^i^^i^^i^^i^'' Such 
flexibility allows for envisioning a vast range of appli- 



cations in optoelectronics (single-electron transistorSji^ 
lasers^SiS£ micro-heaters and micro-refrigerators based 
on thermoelectric effectsSi), life sciences as well as 
in several quantum information processing schemes in the 
solid-state environment i^^iSS*^ 

Almost all QD-based applications rely on (or are influ- 
enced by) electronic correlation effects, which are promi- 
nent in these systems. The dominance of interaction 
in artiflcial atoms and molecules is evident from the 
multitude of strongly correlated few-electron states mea- 
sured or predicted under different regimes: Fermi liquid, 
Wigner molecule (the precursor of Wigner crystal in 2D 
bulk), charge and spin density wave, incompressible state 
reminescent of fractional quantum Hall effect in 2D (for 
reviews see Refs. lllTl . 

It is important to understand where the dominant cor- 
relation effects arise from. In QDs the "external" con- 
flnement potential originates either from band mismatch 
or from the self-consistent fleld (SCF) due to doping 
charges. In both cases, the total field modulation which 
confines a few free carriers is smooth and can often be ap- 
proximated by a parabolic potential in two dimensions. 
Therefore, the single-particle (SP) energy gaps of the free 
carrier states are uniform in a broad energy range. This is 
in contrast with natural atoms, where the electron con- 
finement is provided by the singular nuclear potential. 
Furthermore, while the kinetic energy term scales as rj^, 
rs being the dimensionless parameter measuring the av- 
erage distance between electrons, the Coulomb energy 
scales as . Contrary to natural atoms, in semicon- 
ductor QDs the Coulomb-to-kinetic-energy ratio can be 
rather large (even larger than one order of magnitude), 
the smaller the carrier density the larger the ratio. This 
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causes Coulomb correlation to severely mix many dif- 
ferent SP configurations, or Slater determinants (SDs). 
Therefore, we expect that any successful robust compu- 
tational strategy in QDs, if flexible enough to address 
different correlation regimes, must not rely on trunca- 
tion of the Fock space of SDs obtained by filling a given 
SP basis with N electrons. 

The theoretical understanding of QD electronic states 
in a vastc class of devices is based on the envelope func- 
tion and effective mass modeli2L2& Here, changes in the 
Bloch states, the eigenfunctions of the bulk semiconduc- 
tor, brought about by "external" potentials other than 
the perfect crystal potential, are taken into account by a 
slowly varying (envelope) function which multiplies the 
fast oscillating periodic part of the Bloch states. This 
decoupling of fast and slow Fourier components of the 
wave function is valid provided the modulation of the 
external potential is slow on the scale of the lattice con- 
stant. Then, the theory allows for calculating such enve- 
lope functions from an "effective" Schrodinger equation 
where only the external potentials appear, while the un- 
perturbed crystalline potential enters as a rcnormalizcd 
electron mass, i.e. the "effective" value m* replaces the 
free electron mass. Therefore, SP states can be calculated 
in a straightforward way once compositional and geomet- 
rical parameters are known. This approach was proved to 
be remarkably accurate by spectroscopy experiments for 
weakly confined QDsji2*i& for strongly confined systems, 
such as certain classes of self-assembled QDs, atomistic 
methods might be necessaryi^S 

Several different theoretical approaches have been em- 
ployed in the solution of the few-electron, fully interact- 
ing effective-mass Hamiltonian of the QD systemi^ii^ 
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Here, N is the number of free conduction band elec- 
trons (or free valence band holes) localized in the de- 
vice "active" area, e, m*, and k are respectively the car- 
rier charge, effective mass, and static relative dielectric 
constant of the host semiconductor, r is the position of 
the electron, p is its canonically conjugated momentum, 
A{r) is the vector potential associated with an external 
magnetic field. The effective potential V{r) describes the 
carrier confinement due to the electrostatic interaction 
with the environment. Wc ignore the (usually small) Zee- 
man term coupling spin with magnetic field; this small 
perturbation can always be added a posteriori. 

The eigenvalue problem associated to the Hamilto- 
nian has been tackled mainly via Hartree-Fock (HF) 
method, density functional theory, configuration inter- 
action (CI), and quan tum Monte Carlo (QMC) (for re- 
views sec Ref s . Illl #14lllq) . Each method has its own mer- 
its and drawbacks; mean-field methods, for example, are 



simple and may treat a large number of carriers, a com- 
mon situation in many experimental situations. On the 
other hand, CI and QMC methods allow for the treat- 
ment of Coulomb correlation with arbitrary numerical 
precision, and, therefore, represent the natural choice for 
studies interested in strongly correlated (low density or 
high magnetic field) regimes which arc nowadays attain- 
able in very high quality samples. In addition, CI is a 
straightforward approach with respect to QMC, and gives 
access to both ground and lowest excited states at the 
same time and with comparable accuracy, which allows 
for the calculation of various response functions to be di- 
rectly compared with experiments. The main limitations 
within the CI approach, when applied to the strongly 
correlated regimes attainable in artificial atoms, arc that 
truncation of the Fock space — one of the standard pro- 
cedures in traditional CI methods of quantum chemistry 
(the so-called CI with singles, doubles, etci'^°i'^-^i^^i^^i^'^i'^^) 
— generally gives large errors or even qualitatively wrong 
results, and that the Fock space dimension grows expo- 
nentially with N . 

Here we present a newly developed CI code2& with per- 
formances comparable to QMC calculations for a broad 
range of electron densities, covering the transition from 
liquid to solid electron phases. The implementation is 
independent of the SP orbital basis and, therefore, of the 
details of the device. We do not assume any truncation of 
the Fock space of SDs generated from the chosen SP basis 
[/m;/-CJ30i31,32,33,34,35 (fCI)]. The large eigenvalue prob- 
lem is reduced to its minimal size by exploiting the sym- 
metry of the effective- mass interacting Hamiltonian 
including square total spin [S'^)& The resulting Hamil- 
tonian, which represents the Coulomb interaction in the 
many-electron basis, is diagonalized via a parallel version 
of the Lanczos algorithm. 

The motivation for developing a new code, instead of 
using CI codes already available in the computational 
quantum chemistry literature, was the difficulty of 
handling such codes, optimized for traditional ah initio 
applications, to treat effective or model problems such 
as a few electrons in artificial atoms. Therefore, our 
code has been developed entirely from scratch and 
specifically designed for quantum dot applications. One 
peculiar feature is that it builds the whole Hamiltonian 
in the many-electron basis and it stores it in computer 
memory (see Sec. ^J. This procedure presently con- 
stitutes a limitation with respect to state-of-the-art CI 

^^g.-,j.^^^j^ j30.31.32.33.34.35.38.39.40.41.42.43.44.45.46.47.48.49.50.51.52.53.54.55.5 

of Quantum Chemistry based on the direct tech- 
nique, which avoids the storage of the Hamiltonian 
(cf. Sec. IVIIII . The implementation of the direct CI 
algorithm is left to future work. 

In order to demonstrate the effectiveness of our code, 
we present a benchmark calculation in the strongly corre- 
lated regime par excellence, namely the crossover region 
between Fermi liquid and Wigner crystallization, and we 
compare results with available QMC data. Since the FCI 
code gives access to both wave functions and energies of 
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the first excited states, we analize the signatures of crys- 
tallization in the excitation spectrum. We also discuss 
the performance of the algorithm in a parallel environ- 
ment. 

The paper is organized as follows: First the general 
code structure is illustrated (Sec. ^J, then possible SP 
bases and related convergence problems are considered 
fSec. llll|l . Section HVI analyzes the code performance in a 
parallel environment. A benchmark physical application 
in demanding correlation regimes is considered in Sec.lVl 
where results arc compared with QMC and CI data taken 
from literature. Section IVII discusses the peculiar fea- 
tures of the excitation spectrum in the strongly corre- 
lated regime. Finally, some general considerations are 
addressed fSec. IVIlll . followed by the concluding section 
IVIIII Appendix IXI reviews the construction of the many- 
electron Hilbert space and related issues, while Appendix 
|B]provides the explicit expression of the 2D Coulomb ma- 
trix elements for the so-called Fock-Darwin SP orbitals. 



II. CODE STRUCTURE 

The interacting Hamiltonian can be written in sec- 
ond quantized form-? on the basis of a complete and or- 
thonormal set of SP orbitals {(l)a{r)}: 

^ = J^^'^'^^^'^'^^b'' + \^Y1 ^<^bcdci„cl^,c^ycda- (3) 

ab(T abed crcr' 

Here (caa) creates (destroys) an electron in the spin- 
orbital 0a('")X(T(s), Xt{s) is the spinor eigenvector of the 
spin z-component with eigenvalue a (a = ±1/2), the la- 
bel a stands for a certain set of orbital quantum numbers, 
Eab and Vabcd are the one- and two-body matrix elements, 
respectively, of the Hamiltonian The latter arc de- 
fined as follows: 

Sab = Jdr<p:{r)Hoir)Mr), (4a) 

Vabcd = Jdrjdr' ra{r)rb{r') 

X r-(t)c{r )Mr)- (4b) 

K\r — r I 

The adopted FCI strategy to find the ground- and 
excited-state energies and wave functions of interact- 
ing electrons, described by the Hamiltonian or (jSJ, 
proceeds in three separate steps: 

step 1 A specific set of A^sp single-particle orbitals, 
{(pair)}, is chosen. A"sp should be large enough to 
guarantee as much completeness of the SP space as 
possible. The parameters Sab and Vabcd for a given 
device [Eqs. are computed numerically only 

once, and then stored on disk. In case (j)a{r) is an 
eigenfunction of the SP Hamiltonian Q the matrix 
Sab acquires a diagonal form, Sab — £a6ab- In a few 



(but important) cases, the SP orbitals are eigcn- 
functions of a particularly simple SP Hamiltonian 
Hq, and analytic expressions for Sab and Vabcd can 
be derived; one example is discussed in Sec. IIII Al 

step 2 The many-electron Hilbert (Fock) space is built 
up. The A^-electron wave function, l^'jv), is ex- 
panded as a linear combination of configurational 
state functions^^iSSiSMi (CSFs), namely prim- 
itive combinations of SDs which are simultane- 
ously eigenvectors of the spatial symmetry group 
of Hamiltonian the square total spin S*^, and 
the projection of the spin along the z-axis Sz'. 

=^a,|^,). (5) 

i 

On such a basis set, the matrix representation of 
the Hamiltonian 10) acquires a block diagonal form, 
where each block is identified by a different irre- 
ducible representation of the spatial and spin sym- 
metry group (the latter is labeled by 5^, Sz). The 
algorithm for building the whole set of CSFs start- 
ing from the SP basis set is described in detail 
in App. A; here we are only concerned with the 
fact that CSFs are linear combination of SDs, 
obtained by filling in, in all possible ways consis- 
tent with symmetry requirements, the 2A"sp spin- 
orbitals with N electrons: 

IV'.) =^6y 1$,), (6) 

j 

the 6y 's being the Clebsch-Gordan coefficients. 
These do not depend on the device, and can be 
calculated once (see App. A) and stored on disk. 

step 3 Hamiltonian matrix elements in the CSF represen- 
tation are calculated. Since several CSFs \4>i) may 
consist of different linear combinations of the same 
SDs \^j), it is convenient to compute and store once 
and for all the non-null matrix elements of Hamil- 
tonian Q between all possible SDs, {^j\H\^j'); 
subsequently, the matrix elements between CSFs, 
\ H\ipi'), can be simply evaluated using the co- 
efficients bij previously stored: 

jj' 

Then, the secular equation 

(n - EnI) |*jv) = 0, (8) 

where X is the identity operator, is solved by means 
of the Lanczos package ARPACK^ designed to 
handle the eigenvalue problem of large size sparse 
matrices. Eventually, eigenenergies En and eigen- 
vectors Qi are found and saved, ready for post- 
processing (wave function analysis, computation of 
charge density, pair correlation functions, response 
functions, etc). 
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FIG. 1: Pictorial synopsis of the DONRODRIGO suite of FCI 
codes. 



step 1-3 are implemented in the suite of programs 
DONRODRIGO)^ named after the characters of Alessandro 
Manzoni's Uterary masterpiece / promessi sposi. Each 
step corresponds to a different high-level routine, pro- 
viding maximum code flexibility. Specifically, step 1 
is an independent program on its own, the INNOMINATO 
code: The idea is that, to describe different QD devices, 
it is sufficient to modify only the form of the confine- 
ment potential, V{r), and field components appearing 
in the SP Hamiltonian Indeed, the confinement po- 
tential and the SP basis set only affect parameters Sat 
and Vabcd- Such parameters are passed, as input data, 
to step 2-3 (cf. Fig.nj which constitute the core of the 
main program, and do not depend on the form of the 
SP Hamiltonian and basis set. step 2-3 are by far the 
most intensive computational parts and are parallelized 
via the MPI interface: the Hamiltonian matrix elements 
(<I>j| 7i |<i>j') are distributed among all processors. Only 
non-null elements are allocated in local memory, and the 
parallel version of ARPACK routine is used for diagonal- 
ization. Matrix element computation is optimized by im- 
plementing a binary representation for CSFs: each CSF 
corresponds to a couple of 64-bit integers, where each bit 
represent the occupation of a SP orbital (App. A2), and 
highly efficient bit-per-bit logical operations are used to 
determine the sign of the matrix element, depending on 
the anti-commuting behavior of fcrmionic operators. An 
example is given in Sec. A3 of the Appendix. 

Figure ^ is the graphical synopsis of the interface 
between codes belonging to the suite DONRODRIGO, in- 
cluding post-processing routines, which actually consi- 
tute a rich collection of codes capable of calculating 
experimentally accessible quantities and response func- 
tions. Presently, we have already implemented wave 
function analysis, computation of charge density and pair 
correlation functionsjSiiiiiSiS dynamical form factor^ 
purely electronic Raman excitation spectrum^Sii^ single- 
electron excitation transport spcctrumjiiS spectral den- 



sity weight (one-particle Green's function) jiiSiS^ re- 
laxation time of excited states via acoustic phonon 
emissiouiSSiSl 



III. SINGLE PARTICLE BASIS 

So far, in applications, among all possible SP basis sets 
{(pail")} we chose the eigenfunctions of the SP Hamilto- 
nian HqIt). In this section we review the two most com- 
mon cases. 



A. Fock-Darwin states 

Two common classes of QD devices are the so-called 
"planar" and "vertical" QDs. Planar QDs are obtained 
starting from a quantum well semiconductor hetcrostruc- 
turc, where conduction electrons arc free to move in the 
plane parallel to the interfaces, therefore realizing a two 
dimensional electron gas. The second step is to deposit 
electrodes on top of the structure, which laterally de- 
plete the electron gas forming a puddle of carriers, i.e., 
the do^. By appropriate engineering the top gates, more 
than one QD can also be formed, coupled by quantum 
mechanical tunneling. Vertical QDs are included in cylin- 
drical mesa pillars built by laterally etching a quantum 
well or a heterojunction&i^ In this case, devices with 
two or more tunnel coupled QDs can be obtained start- 
ing from coupled quantum wells or heterojunctions. 

In the above devices, the in-plane confinement is much 
weaker than the confinement in the underlying quantum 
well or hetcrojunction. It is commonly acceptedi^ that 
the confinement potential for the above devices is accu- 
rately described by 
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where luq is the natural frequency of a 2D harmonic trap 
and V{z) describes the profile of a single or multiple 
quantum well along the growth axis z. The eigenfunc- 
tions of the Hamiltonian with the potential given 
by and possibly a homogeneous static magnetic field 
applied parallel to z, are of the type 



(10) 



namely the motion is decoupled between the x — y plane 
and the z axis. Here, the (^„m's, known as Fock-Darwin 
orbitalsi (see below), are eigenfunctions of the in-plane 
part of the SP Hamiltonian with eigenvalues 
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(11) 



n and m being respectively the radial and azimuthal 
quantum numbers [n = 0,1,2,... m = 0, ±1, ±2, . . .), 
fl = {ujQ+uj'^/4y^'^ , LUc ~ eB/m*c, B being the magnetic 
field, and the Xi(z)'s {i = 1,2, 3, . . .) are eigenfunctions 
of the confinement potential along z. At zero field the 
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FIG. 2: Shell structure of the SP Fock-Darwin energy spec- 
trum, and cissociated degeneracies. Each shell corresponds to 
the energy TujJo{Na-hei\ + 1), where A'shcii — 2n + \m\ is fixed, 
and (n, m) are the radial and azimuthal quantum numbers, 
respectively. The degeneracy of each shell (not taking into 
account the spin) is A''sheii + 1- 



Fock-Darwin energies (|llf) display a characteristic shell 
structure, the orbital degeneracy increasing linearly with 
the shell number (Fig.El). 

The explicit expression of Pock Darwin orbitals is 



TT (n + |to|) 



(12) 



where (p, d) are polar coordinates and lL™'(^) is a Gener- 
alized Laguerre Polynomial. The typical lateral extension 
is given by the characteristic dot radius £ = (?i/m*a;o)^^^, 
^ being the mean square root of p on the Pock-Darwin 
ground state iy9oo- 

In usual applications only the first one or two con- 
fined states along z arc important, since the localiza- 
tion in the growth direction is much stronger than in the 
p^g^j^g^9ji3ji4ji^ and it is often possible to assume a mirror 
symmetry with respect to the plane crossing the origin^ 
V{z) = V{—z). In the above simple model the spa- 
tial symmetry group of the device, Dock, is abelian and 
therefore SDs (CSPs) transforming according to its irre- 
ducible representations may be straightforwardly builtiS^ 
The presence of a vertical magnetic field makes the sys- 
tem chiral and reduces the symmetry to Cook- 

The most obvious advantage for choosing as SP ba- 
sis orbitals the eigenfunctions of the SP Hamiltonian it- 
self is that they represent a natural and simple starting 
point with regards to the physics of the problem, allow- 
ing for both fast convergence in SP space and easy sym- 
mctrization of many-electron wave functions. Besides, in 
the two dimensional case Coulomb matrix elements are 
known analytically (cf. App. B), while in the three di- 



mensional case [Eq. ^] one can limit the effort of a six 
dimensional integral evalutation to only a three dimen- 
sional numerical integration, two variables being along 
the z-axis and one along the radial relative coordinate, 
respectively, after transformation to cylindrical coordi- 
nates in the in-planc centcr-of-mass frame. Therefore, 
the evalutation of the Coulomb integrals (|4bp is easily 
feasible — the serial version of the INNOMINATO routine 
contributes a negligible part to the total computation 
time — and makes the usage of Pock-Darwin orbitals 
in the QD context competitive with the implementation, 
e.g., of a gaussian basisiS2i2£ 



B. Numerical single particle states 

Exploiting spatial symmetry and nearly parabolic 
confincmcnte is not always possible. For example, 
self-assembled QDs grown by the Stranski-Krastanov 
mechanism in strained materials, like in several III-V 
materials^i, grow in highly non symmetric shapes. An- 
other case where symmetry is lost is by application of 
a magnetic field with arbitrary directioniSiZSiS In both 
cases, the no analytic solutions are available. In these un- 
favorable cases we solve in a completely numerical way 
the SP Schrodinger equation and obtain Coulomb matrix 
elements via Pourier transformation of Il4bp in the recip- 
rocal space (for more details see Ref. [t^)- Note that the 
Coulomb integrals in l|4b|l become complex as far as the 
field is tilted with respect to the z axis. In this general 
case the algorithm is computationally demanding already 
at the SP level, requiring extensive code parallelization. 



C. A test case — Convergence issues 

In order to assess quantitatively both completeness and 
convergence rates for the Pock-Darwin SP basis, we set 
up a test case, which we will systematically refer to in the 
remainder of the paper. We therefore consider N elec- 
trons in a 2D harmonic trap [Eq. Q with V{z) = 0] as 
the density is progressively reduced, namely the lateral 
dimension i increases as N is kept fixed. This test case 
is a significant benchmark since: (i) path integral QMC 
calculations are available providing "exact" data 
concerning energies and wave functions of the ground 
state; (ii) as far as £ increases, the dot goes progressively 
into a strong correlation regime, which is more and more 
demanding with respect to the size of the many-electron 
space needed; (iii) the model is simple and the parame- 
ters of Eqs. Q may be computed analytically. 

We measure the strength of correlation by means of the 
dimensionless parameter— A = i/a^ [og = fi^ k / {m* e^) 
is the effective Bohr radius of the dot], which is the QD 
analog to the density parameter Ts in extended systems. 
As a rough indication, consider that for A w 2 or lower the 
electronic ground state is liquid, while well above A = 8 
electrons form a "crystallized" phase, rcminescent of the 
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IV. CODE PARALLELIZATION AND 
PERFORMANCES 
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FIG. 3: Relative error of the ground state energy for the 
test case defined in Sec. 1111 CI as a function of the size of 
the SP basis, for different values of the number of electrons, 
A^, and the dimensionless interaction strength. A, for a two 
dimension al harm onic trap. The reference values are taken 
from Refs. l75l76l 



Wigner crystal in the bulk. In the latter phase electrons 
are localized in space and arrange themselves in a geo- 
metrically ordered configuration such that electrostatic 
repulsion is minimized^ 

In Fig. Owe study the completeness of the SP basis, 
plotting the relative error of the ground state energy, as a 
function of the basis size, in different correlation regimes. 
The "exact" reference values are taken to be the QMC 
data^^i^ and the calculated values correspond to adding 
successive Fock-Darwin orbital shells of increasing energy 
to the SP basis (for a total amount of 3, 6, 10, 15, 21, 28 
orbitals, respectively). We see that, in the liquid phase 
(A = 2), a basis made of 15 orbitals is sufficient to guar- 
antee a precision of one part over one hundred. The er- 
ror, however, depends sensitively on A^: the smaller the 
electron number, the higher the precision. The reason 
is that, since the confinement potential is soft, electrons 
tend to move towards the outer QD region as N increases, 
and therefore higher energy orbital shells are needed to 
build wave functions with larger lateral extension. A ba- 
sis made of 21 orbitals, however, is enough to obtain a 
similar precision for both two- and five-electron ground 
states (around 0.1%). Figure |21 also shows that, as we 
move towards the high correlation regime, i.e., as A in- 
creases, the error rapidly increases if the SP basis size is 
kept fixed. For example, when A = 8 the error for the 
N = 5 ground state using 21 orbitals is about 2%, an 
order-of-magnitude larger than for the liquid regime at 
A = 2. Therefore, A turns out to be a crucial parameter 
for the convergence of the SP basis. 



Since the main program of DONRODRIGO (here referred 
to as DONRODRIGO itself) is the most computationally in- 
tensive part of the suite (cf. Sec. EJ, it has been par- 
allelized to take advantage of parallel architectures. In 
this section the parallelization strategies and the results 
of parallel benchmark are discussed. 

As anticipated in Sec.^ the code has been parallelized 
using the message passing paradigm based on MPI inter- 
face. This choice has been guided both by the possibility 
of having a portable code, that could run from small de- 
partmental clusters to large supercomputers, and, most 
of all, by the possibility of using the PARPACK library, 
namely the parallel version of the ARPACK librarjiS 
used in the scalar code to solve the eigenproblem, which 
in turn uses MPI. The package is designed to compute 
a few eigenvalues and eigenvectors of a general square 
matrix, and it is most appropriate for sparse matrices 
(see the ARPACK user guidoS2>). The algorithm reduces 
to a Lanczos process (the Implicitly Restarted Lanczos 
Method) in the present case of a symmetric matrix, and 
only requires the action of the matrix on a vector. For 
a comparative analysis of the algorithm and its perfor- 
mances with regards to other methods see Ref. [771 

Specifically, DONRODRIGO has been parallelized as fol- 
lows: In the first part of the code the matrix elements 
of the Hamiltonian {tj;i\H\ipi') to be computed are dis- 
tributed to the available processes, which then compute 
all the elements in parallel. Once the distributed Hamil- 
tonian matrix has been built, the eigenvalues and eigen- 
vectors are computed using the PARPACK library. Note 
that in the computation of the Hamiltonian no commu- 
nication is needed between processors. The PARPACK 
library manages the communications between processors 
by its own, and the calling code DONRODRIGO has only 
to handle the distribution of data as required by the 
PARPACK itself. Together with the data distribution, 
PARPACK requires that the user writes a few parallel 
support subroutines. Specifically, the computation of the 
product between the Hamiltonian matrix and a vector is 
needed. Within this scheme both data and computations 
are distributed among all processes, and then the parallel 
version of DONRODRIGO scales as the number of processors, 
in terms of both wallckock time and local memory allo- 
cation. 

To evaluate the performance of the code a benchmark 
has been run on two parallel machines, with very different 
architectures: a Linux cluster (CLX, equipped with IBM 
x335, PIV 3.06 GHz, dual processor nodes and Myri- 
com interconnection) and an IBM AIX parallel super- 
computer (SP4, equipped with p690 nodes, containing 
32 Power4 1.3 GHz processors each, and IBM "Federa- 
tion" high performance switch interconnections). 

The test case chosen for the benchmark is the prob- 
lem of iV 4 and = 5 electrons in a 2D harmonic 
trap (see Sec. IIII Cll . with iVgp = 36. Requirements of 
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c 


subspace 


non-null 


size 







8018 


3338976 


1 


11461 


7090087 


2 


3493 


709953 



TABLE I: Subspace dimensions and non-null matrix elements 
for the M = sector of the A*" = 4 problem with Nsp ~ 36. 
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FIG. 4: Overall execution time of DONRODRIGO on CLX and 
SP4, as a function of the number of processors, for a test case 
with iV = 4, Nsp = 36, M = 0, Nsd = 22972. The size 
and number of non-null matrix elements {xpil Ti.\tpii) of each 
subspace is given in Table |I| 

such test in terms of memory and computations arc far 
from typical nee ds of our most demanding applications 
(see, e.g., Refs. I64l65l and Sec. 01. A first small test 
calculation has been chosen where the effect of commu- 
nications over the computations is maximized. It is the 
case of the M = subspace for = 4, where the initial 
SD space with linear size A'sd — 22972, corresponding 
to Sz — 0, is rearranged giving three CSF subspaccs for 
S = 0,1,2. The sizes and numbers of non-null matrix el- 
ements ("0^1 7i iV'i') are given in Table The job is small 
enough to run on a single processor. Even if commu- 
nications between processes dominate, nevertheless the 
execution time (Fig. ^ displays a good scalability up to 
64 processors on both machines, giving the user a large 
flexibility in the choice of the number of processors. 

The speedup (see Fig. ISJ of the computation of the 
Hamiltonian elements (H) scales almost linearly as ex- 
pected, since there arc no communications but those used 
to collect the data. For 64 processors the speedup shown 
in Fig. [SI is sublinear because the data distribution across 
processors turns out to be unbalanced, since the size of 
the job becomes unrealistically small. The computation 
of the eigenvalues and eigenvectors (L) for this test case is 
much less expensive than the computation of the Hamil- 
tonian (roughly 1:10), and its scalability (especially on 
CLX) is very small, since the execution time is already 
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CPUs 

FIG. 5: Speedup vs number of CPUs (ratio of the execution 
time when running on one processor to the time when using x 
CPUs) on CLX and SP4, for the computation of the Hamilto- 
nian elements (H), including I/O, and for the computation of 
eigenvalues and eigenvectors using PARPACK libraries (L), 
respectively. Same test case as in Fig. |1| 



S subspace size non-null 7i 1?/)^') 

M = 9M = M = 9 M = 
1/2 63077 123418 72197095 1.889862-10^ 
3/2 46062 91896 46136336 1.118364-10* 
5/2 9896 20370 2507738 6212517 

TABLE II: Subspace dimensions and non-null matrix ele- 
ments for the M = and M — 9 sectors of the N = 5 Hilbert 
space with Asp — 36. 



very small on one processor. 

In order to evaluate the parallel performance as the 
computational load of jobs increases, we consider two ad- 
ditional test jobs for N ^ 5, corresponding, in order of 
increasing computational load, to M = 9 and AI = 0, 
respectively. In both cases DONRODRIGO rearranges the 
SD space into three CSF subspaces corresponding to 
S = 1/2,3/2,5/2, respectively. The size of each CSF 
subspace and the number of non-null Hamiltonian ma- 
trix elements is given in Table ITTl In Fig.|^the execution 
time on eight SP4 CPUs for the three jobs mentioned 
above is displayed as a function of the job size. The ex- 
ecution time as a function of the linear size of the SD 
space scales roughly as a power law, with exponent close 
to two. Note that both the computation of the Hamilto- 
nian matrix (H) and the solution of the eigenproblem (L) 
have similar behavior. This is possible because the Lanc- 
zos method is applied, which scales with a lower exponent 
than a direct diagonalization method. The ratio between 
the eigenproblem part (L) and the Hamiltonian part (H) 
slowly decreases as the size of the system increases. This 
implies that the parallelism becomes progressively more 
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FIG. 6: Execution time on 8 SP4 CPUs vs SD space lin- 
ear size for selected jobs (see text and Table HTjl . The three 
lines refer to the overall execution time (Total), the time 
spent to compute Hamiltonian matrix elements (H) , including 
I/O, and that to compute eigenvalues and eigenvectors using 
PARPACK libraries (L). 



effective as the system size increases, since part H dis- 
plays a better speedup than part L when the number of 
processors increases (cf. Fig. EJ. 



V. TWO-DIMENSIONAL QUANTUM DOT IN 
VERY STRONG CORRELATION REGIMES 

In this section we demonstrate the power 
and rehability of DONRDDRIGD in treating the 
few-electron problem in QDs for our test case 
(cf. Sec. ImC)! . There have been many CI calcula- 
tions for few-electron QDs in different geometries and 

conditions .R-3-11-12.1.T14.1R.fi.Tfi4.fifl.fi3.7n.72.78.73.Sn.S1.S2.8.S.84.8R.! 

Here we consider iV electrons in a 2D harmonic trap, as 
the density is progressively diluted, namely going from 
a Fermi liquid behavior (small A) to a Wigner molecule 
regime (large A). Such a case is particularly interesting 
since it is the object of a querelle between different 
theoretical methods (HF^ QMCr^^ C&.) about what 
value of A corresponds to the onset of "crystallization" 
and what the spin polarization of the six-electron ground 
state is (see the introduction of Ref. for a discussion 
of the physics behind this problem). 

A complete analysis of the physics emerging from our 
results is beyond the aim of this paper. In this section 
we just focus on the convergence and accuracy of ground 
state energies, for different and Hilbert space sectors, 
in a wide range of A. A possible source of confusion 
in the above querelle is the fact that different authors 
use different conventions to identify density values, like 



;74^76,95^96 q^. ^^^93J^ qj. other i2i Therefore, it is not al- 
ways clear which energy values, taken from different pa- 
pers, should be compared. Here we consider path integral 
QMC datapiiS^ for 2 < iV < 8 and 2 < A < 10 and CI 
data obtained from a cutoff procedure applied to the SD 



A M S 



CM 



6 

1 



8 

1 



10 

1 



CO 



10 1 





6 


2 



10 


2 



SP set dimension 



21 



28 



36 



3.7338 3.7312 3.7295 
4.1437 4.1431 4.1427 



4.8502 
5.1203 



5.7850 
5.9910 



6.6185 6.6185 6.6185 
6.7873 6.7873 6.7873 



7.3840 
7.5286 



21 



36 



3.1755 8.1671 
3.3244 8.3224 



11.046 11.043 
11.055 11.053 



13.468 13.467 
13.439 13.438 



15.641 15.634 
15.597 15.595 



17.652 17.630 
17.600 17.588 



36 

13.626 
13.771 
13.848 
14.256 



19.035 
19.146 
19.170 
19.359 



23.641 23.598 
23.691 23.650 
23.870 23.805 



32 



36 



27.677 27.675 
27.702 27.699 
27.828 



31.429 
31.441 
31.553 



Rcf. 74.75 Ref. 76 



4.893(7) 
5.118(8) 



8.16(3) 
8.37(1) 



11.05(1) 
11.05(2) 



11.055(8) 
11.050(10) 



13.43(1) 



15.59(1) 



17.60(1) 



13.78(6) 



14.30(5) 



19.15(4) 



19.42(1) 



19.104(6) 
19.34(1) 



23.62(2) 
23.790(12) 



27.72(1) 
27.823(11) 



31.48(2) 
31.538(12) 



TABLE III: Comparison between FCI ground-states energies 

i.foateoofc Mv.QiM^Mmo^i . <w.%8 Wfi*MGfb2a«i . BhlQ, . fe§ afeio8.io9.iio.i] 

taSen from the literature for 2 < A^ < 4 electrons in a m 

harmonic trap. Few-electron states belong to different Hilbert 
space sectors classified according to their total orbital an- 
gular momentum hA4 (z-component) and total square spin 
h'^S{S + 1) (units of 1/2). The parameter A is the ratio be- 
tween the characteristic dot radius £ and the effective Bohr 
radius (see text). In the middle column, we report results cor- 
responding to diferent SP basis set, as indicated in the cells 
with gray background. 



space2^ for = 3, 4 and 2 < A < 20. Note that the 
extreme values of density and N considered here cannot 
be reached simultaneously, since they arc well beyond the 
limits of FCI calculations performed until now for QDs. 
For example, the lower density limit for the A^ = 6 CI 
calculation of Ref. HI was A « 3.5, while Mikhailo\i2^ 
considered a large density range but only up to A^ = 4. 

In Tables IIIII and IIVI we compare our FCI results for 
ground state energies, for 2 < A < 10 and 2 < A^ < 8 
and different values of 5* and M, with QMC data taken 
from Refs. The energies are given in units of 

Tlujq. Progressively larger sets of SP orbitals are consid- 
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\ M S 


SP set dimension 


Ref . 7^^^ Ref . ^6 




2 11 

2 3 
1 3 
5 


21 28 36 

20.36 20.34 20.33 
20.64 20.62 20.61 
20.92 20.90 20.90 
21.15 21.13 21.13 


20.30(8) 
20.71(8) 

21.29(6) 




4 11 

2 3 

5 


29.01 28.96 28.94 
29.21 29.11 29.10 
29.44 29.31 29.30 


29.09(6) 29.01(2) 
29.15(6) 29.12(2) 
29.22(7) 29.33(2) 


6 11 

2 3 
5 


36.22 
36.34 
36.42 


36.26(4) 
36.35(4) 
36.44(3) 


8 11 

2 3 
5 


28 32 36 

42.98 42.80 42.77 
43.04 42.91 42.86 
43.01 42.93 42.88 


42.77(4) 
42.82(2) 
42.86(4) 




10 1 1 

1 3 

2 3 
5 


48.84 
48.91 
48.93 
48.91 


48.76(2) 
48.78(3) 

48.79(2) 




2 
1 2 
4 
6 


21 28 36 

28.03 27.98 
28.36 28.30 
28.54 28.48 
28.98 28.93 




CO 


4 
1 2 

4 
6 


40.74 40.54 40.45 
40.96 

41.06 40.71 40.66 
41.29 40.88 40.85 


40.53(1) 
40.62(2) 
40.69(2) 
40.83(4) 


6 

4 
6 


51.35 51.02 
51.38 51.15 
51.46 51.25 




8 
1 2 
4 
6 


26 28 36 

61.44 61.35 60.64 
61.53 61.37 60.71 
61.43 61.32 60.73 
61.47 61.38 60.80 


60.37(2) 
60.42(2) 




10 

4 
6 


69.74 
69.81 
69.86 






4 3 7 

5 

1 3 
1 

2 1 


15 21 

57.02 55.20 
56.53 54.93 
56.55 54.78 
56.60 54.69 
54.68 


53.93(5) 
53.80(2) 

53.71(2) 


00 

II 


2 2 


4 
6 
8 


15 21 

47.14 
47.26 
47.58 
48.19 
48.48 


46.5(2) 

46.9(3) 
47.4(3) 
48.3(2) 


4 2 

4 
6 

8 

1 
1 2 


73.54 70.48 
73.63 70.57 
73.72 70.65 
74.17 70.73 
74.15 70.83 
73.44 
73.44 


68.44(1) 
68.52(2) 

68.3(2) 
68.5(2) 
69.2(1) 

68.44(1) 



TABLE IV: Same comparison as in Table [TTTI for 5 < iV < 8 
electrons. 



ered. Note that, contrary to available QMC energies, 
which only refer to different values of S, here in addi- 
tion we are able to assign values of M. QMC data for 
a given value of S are compared with FCI data with M 
corresponding to the lowest energy for the same value 
of S. For 2 < iV < 5 we find an almost perfect agree- 
ment between FCI and QMC data in the full range of A; 
in several cases FCI results are slightly lower in energy 
than corresponding QMC data. Since the FCI calcula- 




dimensionless interaction strengtii X 



FIG. 7; Weight of the CSF with the largest probability in 
the linear expansion of the interacting ground state, in a two 
dimensional harmonic trap, as a function of the dimensionless 
interaction strength A, defined in the main text, for different 
values of the number of electrons, A^. Note that, for N = 3, 
since there is a crossing in energy between singlet and triplet 
lowest-energy states in the range 4 < A < 6, both states are 
depicted. 



tion takes into account all possible SDs that can be made 
from a given SP orbital set, the accuracy of the calcula- 
tion solely depends on the size of the SP basis. In this 
regard, we note that a SP basis with A^sp = 21 is large 
enough to achieve full convergence of results, at least for 
A^ < 4. For A^ = 5, a good compromise is A^sp = 28 
fTable lTvji . For A^ > 6, the larger A, the larger the re- 
quired A'sp to achieve convergence. For A^ = 6 we were 
able to obtain well converged results in the full A range 
only using A'sp = 36. The agreement is excellent up to 
A « 8 (the slight discrepancy between FCI and QMC 
data is five parts per thousand for A = 8). For N > 7 
we were only able to use A'sp = 21, and of course the 
FCI-QMC discrepancy gets worse. However, even in this 
case, the agreement between FCI and QMC data is at 
least qualitatively correct. In fact, FCI and QMC agree 
on the value of 5* for the absolute ground state. Indeed, 
a wrong prediction for S would be a typical indication of 
bad convergence. 

The above discussion shows that: (i) FCI gives results 
with accuracy comparable to that of QMC (and in addi- 
tion provides the wave function of both ground and ex- 
cited states) (ii) FCI calculation becomes more demand- 
ing as A (and A^, of course) increases. 

To gain a deeper insight into point (ii), we analyzed the 
ground state wave function for different values of A and 
A^, focusing on the CSF with the largest weight in the 
FCI linear expansion (O of the interacting wave function. 
Figure [7| shows that the weight of the dominant CSF 
dramatically decreases as either A or A^ increase. For 
example, at A = 2, the dominant CSF weight changes 
from w 80% for A^ = 2 to « 20% for A^ = 5. The weight 
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loss is even more dramatic as A increases: e.g., for iV > 4 
the decrease is larger than one order of magnitude as 
A goes from 2 to 10. This trend clearly illustrates the 
evolution between two limiting cases: as A — *■ 0, the dot 
turns into a non-interacting system, and just one CSF is 
the exact eigenstate of the ideal Fermi gas, with weight 
equal to 100%. What CSF is the ground state is dictated 
by the Aufbau principle, according to which the lowest- 
energy SP orbitals are filled in with N electrons. As 
X —^ oo, the system evolves into the limit of complete 
crystallization: electrons are perfectly localized in space, 
therefore the SP basis, and consequently the CSF basis, 
turn out to be extremely ineffective in representing the 
Wigner molecule, unless one uses a special SP basis of 
localized gaussians, or similar kinds of orbitals. 

The loss of weight shown in Fig. [7| depends of course 
on the chosen SP basis. Therefore we expect that, by 
means of a proper unitary transformation of SP orbitals, 
like that provided by a self-consistent HF calculation, 
weights of dominant CSFs will be in general larger than 
those shown in Fig. |7| However, except for a shift of 
weight values, the above discussion will not be qualita- 
tively affected. Note that past experience in HF calcula- 
tions showed that convergence of the HF self-consistent 
cycle is not necessarily achieved at low densitiesiiS^ 

Figure [7| also shows that the triplet ground state, for 
iV = 3, is relatively less affected by interaction than the 
singlet ground state (see also Fig. |HJ), namely less CSFs 
are needed for the triplet than for the singlet in the FCI 
linear expansion |SJ) of the wave function. The reason 
is that in the triplet exchange interaction is sufficient to 
keep electrons far apart as Coulomb repulsion between 
electrons increases; the same effect can be obtained for 
the singlet only in virtue of a deeper correlation hole 
around each electron. Such behavior is obtained by mix- 
ing more CSFs, namely the singlet increases its correla- 
tion with respect to the triplet. 

Not only the dominant CSF loses weight as A — > oo, 
but also more and more CSFs are needed to build up 
the exact wave function. This is illustrated by Fig. |S1 
where we show the sum of weights of "significant" CSFs 
as a function of A and for different values of N. Here 
by "significant CSF" we mean a configuration \ipi) whose 
weight \ai\ is larger than 1% [cf. Eq. ©]. Figure |S| 
clearly demonstrates that, as A increases, the set of sig- 
nificant CSFs is progressively emptied. From such result 
we infer that, in regimes of strong correlation, only the 
full CI approach is a reliable al gori thm. Therefore, cut- 
off procedures like those of Ref. |93, where the full space 
of SDs is filtered and truncated by considering the ki- 
netic energy per single SD, are potentially dangerous in 
strong correlation regimes where they might not be able 
to warrant convergence si^E 

This is illustrated in Table |3 where we compare CI 
results obtained via two different algorithms. The first 
column in Table shows ground state energies, for 
2 < A < 20 and = 3, 4, obtained from FCI DONRODRIGO 
runs. Here A^sp = 55 was considered, and no cutoff pro- 




dimensionless interaction strength X 



FIG. 8: Total weight of CSFs whose probability in the linear 
expansion of the ground state is larger than 1% as a function 
of the dimensionless interaction strength A for different values 
of A'^. The system is the same as in Fig.Q 
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3 


8.1633 
8.3217 


8.1651 
8.3221 


CO 

II 


4 11 

3 


11.0423 
11.0527 


11.0422 
11.0527 


6 11 
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13.4666 
13.4380 


13.4658 
13.4373 


8 11 
3 


15.6344 
15.5948 


15.6334 
15.5938 




10 1 1 

3 


17.6293 
17.5877 


17.6279 
17.5863 


15 1 1 

3 


22.1127 
22.0754 




20 1 1 

3 


26.1184 
26.0863 






2 2 

2 4 


13.6195 
14.2544 


13.6180 
14.2535 


II 


4 2 

2 4 


19.0335 
19.3578 


19.0323 
19.3565 


6 2 

2 4 


23.5975 
23.8041 


23.5958 
23.8025 


8 2 

2 4 


27.6715 
27.8222 


27.6696 
27.8203 




10 2 

2 4 


31.4148 
31.5352 


31.4120 
31.5323 


15 2 

2 4 


39.8197 
39.9054 


39.8163 
39.8970 


20 2 

2 4 


47.3443 
47.4153 


47.4002 
47.4013 



TABLE V: Comparison between different CI calculations for 
— 3,4. FCI ground-state energies (units of hujo) from 
DONRODRIGO were obtained using a SP basis set corresponding 
to the first ten lowest-energy Fock-Darwin shells (55 orbitals) 
and no truncation of the many-electron space. Results taken 
from R.efs. IqrIIq^ were obtained considering a larger number of 
SP orbitals but with a cutoff procedure on the kinetic energy 
per single SD. S is given in units of 1/2. 
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cedure was implemen ted. T he second column collects 
data taken from Refs. I95l96i In those papers Ngp > 55 
but only certain SDs are retained, namely those whose 
kinetic energy per single SD is smaller than a certain 
cutoff. For example, in Ref. [o^ the largest size of the 
SD space with TV = 4, S* = 1, M = (5 = 2, A/ = 2) 
considered was 24348 (8721) while in our calculation was 
67225 (15659). The comparison between the two sets of 
data gives excellent agreement, since results systemati- 
cally have at least four significant digits identical (ener- 
gies in the second column tend to have a slightly lower 
energy due to the larger A^sp which could be used for 
such a small number of electrons). However, there is one 
important exception, namely at very low density (A = 20, 
TV = 4, 5 = 1, M = 0), where our result is 47.34 against 
47.40. This can be explained by the loss of angular cor- 
relation once one truncates the SD space as it is done in 
Ref.ii 



VI. NORMAL MODES OF THE WIGNER 
MOLECULE 
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FIG. 9: Excitation energies of a six-electron quantum dot in 
the Wigner regime (A = 8) for different values of the total 
orbital angular momentum, M, and total spin, S. Energies 
are in units of huo. 



In order to illustrate the wealth of information pro- 
vided by the FCI method with respect to QMC tech- 
niques, in this section we give an example of excita- 
tion spectrum calculation. In fact, the access to excited 
states is a difficult task for QMC calculations, which rou- 
tinely focus on ground-state properties'^ Nevertheless, 
the knowledge of excited-state energies and wave func- 
tions is of the utmost interest from both the theoretical 
and experimental points of view. Indeed, several useful 
dynamical response functions can be computed from ex- 
cited states. Moreover, many experimental techniques 
allow nowadays to access QD excitation spectrum, such 
as single-electron non-linear tunneling experiments^^*i^ 
far-infrared spctroscopyji inelastic light scatteringii^ In 
addition, the theoretical analysis of the energy spectrum 
gives a precious insight into QD physics. 

Figure^ldisplays our FCI results for the low-energy re- 
gion of the excitation spectrum of a six-electron quantum 
dot in the Wigner regime (A = 8), for different values of 
the total orbital angular momentum and spin multiplic- 
ities. The plot of the lowest energies as a function of 
the angular momentum is known as yrast line in nuclear 
physics^ and provides several hints on the crystallized 
electron phase. 

We see from Fig. |31 that the ground state is the non- 
degenerate spin singlet state with T\/ = 0, which is 
found to be the lowest-energy state in the whole range 
< A < 10. The absence of any level crossing as the den- 
sity is progressively diluted (i.e. A increases) implies that 
the crystallization process evolves in a continuos manner, 
consistently with the finite-size character of the system. 
Nevertheless, several features of Fig. O demonstrate the 
formation of a Wigner molecule in the dot. 

First, the three possible spin multiplets other than 
5 = 0; namely 5* = 1,2,3, lie very close in energy to 



the ground state. In general, at the lowest energies for 
any given M, several spin multiplets in Fig. [^appear as 
almost degenerate. 

The overall behavior is well explained by invoking elec- 
tron crystallization. In the Wigner limit, the Hamilto- 
nian of the system turns into a classical quantity, since 
the kinetic energy term in it may be neglected with re- 
spect to the Coulomb term. Therefore, only commuting 
operators (the electron positions) appear in the Hamil- 
tonian. In this regime the spin, which has no classical 
counterpart, becomes irrelevant spin-dependent ener- 
gies show a tendency to degeneracy. This can also be 
understood in the following way: if electrons sit at some 
lattice sites with unsubstatial overlap of their localized 
wave functions, then the total energy must not depend 
on the relative orientation of neighboring spins. 

Note also that the proximity of the four possible spin 
multiplets at the lowest energies occurs for both T\/ = 
and for M = 5. Such a period of five units on the M axis 
identifies a magic number^^ii^ whose origin is brought 
about by the internal spatial symmetry of the interacting 
wave functiouiiS^ In fact, when electrons form a stable 
Wigner molecule, they arrange themselves into a five-fold 
symmetry configuration, where charges are localized at 
the corners of a regular pentagon plus one electron at the 
centeriSiiSiiaa 

A second distinctive signature of crystallization is the 
appearance of rotational hands^^ which we identify in 
Fig. 1^1 as those bunches of levels, separated by energy 
gaps of about 0.2?ia;o, that increase monotonically as M 
increases. We are able to distinguish in Fig. |51 at least 
three bands, where each band is composed of different 
spin multiplets. Such bands arc called "rotational" since 
they can be identified with the quantized levels E^ot{M) 
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of a rigid two-dimensional top, given by the formula 




where / is the moment of inertia of the topii2i These 
excitations may be thought of as the "normal modes" of 
the Wigner molecule rotating as a whole in the xy plane 
around the vertical symmetry axis parallel to z. 

VII. DISCUSSION 



respect to previous CI codes, DONRDDRIGO is powerful 
enough to treat demanding correlation regimes with the 
same degree of accuracy as QMC method. In addition, 
the code provides full information about ground and ex- 
cited states, its only limitation being the number of elec- 
trons. Its flexible structure allows for inclusion of dif- 
ferent device geometries and external fields, plus easy 
post-processing code development. 

DONRDDRIGO is available, upon request and 
on the basis of scientific collaboration. at 
^http : / / www ■ s3 ■ inf m . it/ donrodr igo 



In the previous sections we demonstrated that 
DONRDDRIGO is a flexible and well performing code, suit- 
able to treat strongly-correlated few-electron problems in 
quantum dots to the same degree of accuracy as QMC 
calculations. With respect to the latter, the FCI method 
also provides full access to excited states. 

A major achievement of DONRDDRIGO is the indepen- 
dence of core routines from both Hamiltonian and SP 
basis. Differently from other CI codes, the SP basis is 
not optimized by a SCF calculation, before the actual 
FCI calculation. Of course the FCI results are indepen- 
dent from the SP basis, but the question arises whether 
a proper SP basis could improve the efficiency of the FCI 
calculation, at the time cost of a preliminary SCF cal- 
culation. However, the usage of non self-consistent SP 
orbitals does not seem to be a serious drawback, since 
DONRODRIGO was designed to treat a wide range of differ- 
ent correlation regimes, especially those for which SCF 
calculations are not expected to converge. An example 
of the latter circumstance is the case of very large dots 
(A large) . Of course in such regime a large SP basis and 
a full CI diagonahzation are both necessary requirements 
(see discussion in Secs. lIIII andlVll. Another major advan- 
tage is the code simple and transparent structure, which 
allowed for easy parallelization and straightforward cod- 
ing of post-processing routines. 

Although very efficient also in highly correlated 
regimes and with very good scaling properties in parallel 
architectures, our code still lacks some smartnesses 
of pre-existing quantum chemistry CI codes. One 
major drawback is the requirement of large amounts 
of memory, since, for each Hilbert space sector, all 
matrix elements {^j\H\^j') are stored (cf. Sec. ITl|l . 
An alternative strategy is given by so-called direct 

30.31.32.33.34.35.38.39.40.41.42.43.44.45.46.47.48.49.50.51.52. 

where required matrix elements are computed on the 
fly during each step of recursive diagonahzation of 
the many-electron Hamiltonian. The latter approach 
constitutes an interesting direction for future work, 
together with possible interfaces with other CI codes. 



VIII. CONCLUSION 

We presented a new high performance full CI code 
specifically designed for quantum dot applications. With 
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APPENDIX A: CONFIGURATIONAL STATE 
FUNCTIONS 

In this Appendix we illustrate how CSFs are built. The 
method is straightforward, namely one writes the op- 
erator in a suitable matrix form and diagonalizes it nu- 
merically, storing the eigenvectors (the Clebsch-Gordan 
coefficients bij) once and for all. This approach is men- 
tioned for example in Ref. |63, together with other possi- 
ble strategies. However, it turns out that the approach 
is non standard for actual CI code implementations^^, 
and we summarize here the algorithm for completeness. 
By means of simple examples we show how it is imple- 
mented in DONRODRIGO, following Slater^ In this section 
we distinguish operators from eigenvalues by means of the 
"hat" symbol. 



53.54.55.56.57.58 



1. Building the configurational state function space 

Let us consider a SP basis made of iVsp orbitals. The 
corresponding SDs are built occupying such orbitals with 
N electrons in all possible ways, taking the two possible 
spin orientations into account. A^sd = C^^) is the total 
number of SDs so obtained; each SD is trivially asso- 
ciated with an eigenvalue of Sz- Note that many SDs 
may be associated to the same value of Sz- The SD set 
!<!>;) = i = 1,2,..., A^sD constitutes a basis for both 
operators Sz and Ti.. Therefore, the energy eigenvalues 
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can be obtained by diagonalization of 



{hSzj\li\i,Sz 



hi = 1, 



,iVsD. 



(Al) 



If Szi 7^ Szj the matrix element (|A1|) is zero, i.e., TL is 
block diagonal in the different Sz sectors. Since both 
operators Sz and commute with the Hamiltonian, as 
well as with each other, it is also possible to build a Fock 
space whose basis vectors are simultaneously eigenstatcs 
of operators 5'^ and 5^. In this new basis, 7i has a block 
diagonal representation, labelled by {Sz^ S), with dimen- 
sions of the blocks smaller than in the previous case. Be- 
low we show how one can find simultaneous eigenstates 
of both S*^ and 5*2. 

The search of the eigenstatcs of 5*^ starting from those 
of Sz can be traced to that of a unitary transformation 

\^)^\^), (A2) 

where \i) are eigenstates of 5*2 and \i) of Sz and S^ simul- 
taneously. It is possible to represent the vectors of the 
second basis on those of the first through a suitable linear 
combination with (Clcbsch-Gordan) coefficients 6^-, 



_ NsD 

under the unitarincss condition 



\J), 



(A3) 



Nsu 

E 



(A4) 



The problem of finding the Clcbsch-Gordan coefficients 
can be greatly reduced by considering only the spins of 
singly occupied orbitals, which, in general, are much less 
than the total number of electrons. This is based on the 
following property: the SDs with only empty or doubly 
occupied orbitals are eigenstates of 5^ with 5 = 0. To 
see this, consider the identity {h = 1) 



S' 



{Sx + iSy){Sx - iSy) + 5^ 



(A5) 



which does not depend on N. Furthermore, let |t) = 
c||0) and ||) = c||0) denote the spinors related to a 
generic orbital. Then the following relations hold: 



{Sx-^Sy)\^) = \i), 
iSx+^Sy)\i) = \^), 

{Sx ~iSy) \i)^0, 

{Sx + tSy)\i;) = 0, 
^.IT> = i/2|T), 

Sz\i) = -l/2\i). 

The operator {S^ — iSy) is a step-down operator: when 
acting on a given SD. the result is a sum of SDs each 



having one of the original up spins turned down, 
example. 



For 



(5,-*5,)|TTi> = UTi> + ITii), (A6) 

where |tTi) = ^a'\'^\-\^\i\^) ■ Analogously, the operator 
(Sx + iSy) is a step-up operator, turning [ spins into t 
spins. The above operators in second quantized form, for 
just one orbital, are: 



{Sx + iSy) = 5^ 



{Sx - iSy) = 5 = c|C| 
5. = 1/2 ' ' 



We now consider two electrons on the same orbital. The 
corresponding SD, c|c| |0), is an eigenstate of Sz with 

Sz = To see that it is also an eigenstate of S^ with 
5 = it is sufficient to apply (|A5|) using the definition of 
S^ and One verifies immediately that 

5+5-44 |0) = c\cic\c^c\c\ |0) = 0. (A7) 

In the general case, N electrons are distributed over A^gp 
orbitals, so that TV < 2A^sp- The operator S is 



S = S(l) + S(2) 



S(A^sp), 



(A8) 



where S(a) = [Sx{a)^ Sy{a)^ Sz{a)\ acts on the a-th or- 
bital. Therefore, we can write 



S' = 



Nsp 
a, 6=1 



5^ Sz 



(A9) 



Again, we ask whether a generic SD, which is eigenstate 
of Sz with some orbitals a doubly occupied or otherwise 
empty, is also eigenstate of 5^ . The answer is affirmative, 
since, in the term 5+5- explicitcd in Eq. IjAQp . when 
first 5- acts on the SD, it destroys a spin f in the full 
a-th orbital and creates a spin | on the same orbital; 
the total contribution of 5+5~, however, is zero, since 
the a-th orbital is already occupied by a spin | electron. 
Therefore, it is proved that the SD is eigenstate of 5^ 
with 5 = 0. 



2. Examples 

Let us consider the case of SDs having some singly 
occupied orbitals. For definiteness, let us focus on the 
example A^ = 4 and A'sp = 4 (a = 0, ... ,3), with two 
electrons in orbital 0, zero in 2. and one in each orbital 
1 and 3, with opposite spin: 



|0) 



(AlO) 



The operator Sz, explicitly, is 



1 ^ 

2E(4tM-4c,x). (All) 



a=0 
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Evidently, Sz\ii) = 0\ii)- The goal is to build linear 
combinations of SDs, including \ii), that are eigenstates 
of S^. To this aim one needs to identify all the SDs which 
couple to |zi) through the operator and represent S*^ 
in this set. In the present case, it is sufficient to consider 
only SDs with Sz = 0. It is possible to establish a priori 
what these SDs arc since they differ from \ii) only for 
the spin orientation of the electrons on the same singly 
occupied orbitals (1 and 3 in the present example), while 
the doubly occupied or empty orbitals (0 and 2, respec- 
tively) remain unchanged. This is evident from Eq. ljA9|l . 
using similar arguments as above; for a formal proof see, 
e.g., Ref. I35I This way, starting from \ii), we build the 
equivalence class S (of which is the representative) 
of all SDs having the same doubly occupied and empty 
orbitals, and differing only in the spin configurations of 
same singly occupied orbitals, with the same value of Sz- 
In our example, the only other SD belonging to S is 



N2) — cLcLcJi cL |0) 



(A12) 



The idea is that we only need to find the Clebsch- 
Gordan coefhcients of the SDs belonging to the same 
equivalence class S. In this way wc reduce the problem 
to the vectorial sum of n angular momenta of value 1/2, 
where n is the number of unpaired electrons, and usually 
n ^ N. We approach this problem by first writing down 
as a matrix, and then diagonalizing it numerically: 
the eigenvectors are the requested Clebsch-Gordan coef- 
ficients. In the case considered above, the two elements 
belonging to S are 



Ti 



it 



(A13) 



where according to the above prescription we only re- 
fer to the spins of singly occupied orbitals, 1 and 3, as 
indicated. On this basis, the matrix form of S'^, using 
Eq. lUnil, is 



(A14) 



whose eigenvalues correspond to the singlet {S — 0) and 
triplet {S = 1) states. The corresponding eigenvectors 
provide us the desired Clebsch-Gordan coefficients 6^, 
i,j = 1,2 [cf. Eq. ljA3|l ]. Specifically, the eigenstates of 
obtained from the SDs belonging to the equivalence 
class S are: 




1 



1 



V2 



1 



V2 



1 



— \^,)- — \^^) (5 = 0). 



^/2 



V2 



(A15) 



(A16) 



Note that, for the triplet state S - 
CSFs degenerate in energy with Sz 



1, there are other 
= ±1, namely |tT) 



and III) (incidentally, the CSFs with S = ±n are always 
single SDs). Such CSFs are redundant in applications 
and are ignored by DDNRDDRIGO. In general, the subspace 
diagonalization with Sz ~ (or Sz = 1/2 if n is odd) 
provides all spin eigenvalues allowed by symmetry. Such 
sectors are the only ones considered by the code. 

As another example, let us consider the case with 
A''sp = 15 (the SP orbital index runs from to 14) and 
ri = 3 singly occupied orbitals {N can be, of course, much 
larger). Let us now focus on a given SD, e.g.. 



|0) 



{N = 7). The above SD corresponds to a state with 
two doubly occupied orbitals (0 and 7), two singly oc- 
cupied orbitals with spin up (3 and the 5), one singly 
occupied orbital with spin down (14). In DONRODRIGO a 
SD is uniquely identified by a couple of 8-byte integers in 
binary representation. The first (second) integer labels 
the spin-up (-down) orbitals: each of the 64 bits repre- 
sents a SP orbital; if a bit is set to 1 then the orbital 
is occupied by one electron. Therefore, the above SD is 
coded as 



12 3 4 



6 7 8 D 10 11 12 13 14 



100101010000000 
1 0000 1 00 1 



(A17) 



The corresponding equivalence class S uniquely associ- 
ated to (|A17|1 is coded as 



12 3 4 



D 10 11 12 13 14 



00010100000000 1 
1 0000 1 00 



(A18) 



where the first integer labels the singly occupied orbitals: 
if a bit is set to 1 then the orbital is occupied by one 
electron. The second integer represents doubly occupied 
orbitals, using the same convention. 

Independently from the number of doubly occupied or- 
bitals, it will be sufficient to solve the equivalent Clebsch- 
Gordan problem for the following three determinants. 



ITTI) UTT) int), 



(A19) 



where the first ket, |tTI)i represents (|A17|I . Note that for 
all kets Sz = 1/2. Therefore, we need to diagonalize S^ 
— a 3 X 3 matrix — as we did with (|A14|) . One obtains 
3 eigenvalues, 1/2,1/2,3/2, and the associated Clebsch- 
Gordan eigenvectors. Eventually, the equivalence class 
(IA18|I is associated with an additional index labelling 
the pertinent Clebsch-Gordan eigenspace (in this case, 
either the quadruplet or one of the two doublets). These 
two entities — the equivalence class (|A18ll and the latter 
multiplet index — identify the CSFs. 

More generally, a A^-electron SD, with n' electrons 
in doubly occupied orbitals and n in singly occupied or- 
bitals {n + n' = iV), is uniquely associated with an equiv- 
alence class Si. The class Si is the same for all and only 
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M S = S^l S = 2 S = 3 NsD S^ 

661300 1131738 568896 97976 2459910 190380 

1 656476 1123952 564697 97221 2442346 189000 

2 643242 1100391 552661 95043 2391337 185624 

3 621112 1062496 532877 91493 2307978 179728 

4 591897 1011245 506545 86741 2196428 172093 

5 555754 949079 474285 80960 2060078 162429 

6 514945 877805 437690 74407 1904847 151662 

TABLE VI: Partitioning of the SD space into CSF subspaces 
for selected values of the total orbital angular momentum M. 
We consider A'' = 6 and a SP basis made of Fock-Darwin 
orbitals corresponding to the first 8 shells {Nsp = 36). S is 
the total spin, Nsd is the dimension of the SD space, and Si is 
the number of equivalence classes from which CSFs are built. 



those SDs which differ from \i) for the spin orientation of 
the n electrons, keeping 5*2 constant. DONRODRIGO scans 
the whole SD space exhaustively (i = 1, . . . , A^sd), 

rearranging it as a set of different equivalence classes Si . 
Within each class iS^, the problem is now equivalent to 
the vectorial sum of n spins, the maximum possible value 
of n (S) being A^ {N/2). For large values of n, where the 
solution of matrices like (|A14p is not possible by sim- 
ple analytical methods, is diagonalized by a standard 
numerical routine. In such a way Clebsch-Gordan coef- 
ficients up to n = 15 (corresponding to a 6435 x 6435 
matrix) were obtained and stored. By means of associ- 
ating each class Si with all possible multiplet indices, all 
CSFs are built. 

As an example taken from an actual 
application^SMSiiSS we consider the reduction of the 
many-electron Hilbert space in our test case (Sec. IIIC) 
of a 2D harmonic trap with A^ = 6. Truncating the 
SP basis to the first 8 shells of Fock-Darwin orbitals 
(A^gp = 36), we first reduce the SD space in blocks 
corresponding to different values of the total orbital 
angular momentum, M (see Table IVI|I . Focusing on 
M = 0, the dimension of the SD space corresponding 
to 5*2 = is A^sD = 2459910. From such starting 
point DONRODRIGO extracts 190380 equivalence classes, 
allowing to build a CSF space divided in singlet {S ~ 0, 
dimension 661300), triplet (5* = 1, dimension 1131738), 
quintuplet (S = 2, dimension 568896), and septuplet 
(5 = 3 dimension 97976) sectors, respectively. The sum 
of sector dimensions is equal to A'sd- Table IVII is a 
larger list of CSF subspaces partitioning the initial SD 
space for different values of M. 



3. Computing matrix elements 

Since any state may be written as a linear combina- 
tion of CSFs I'i/'i), any matrix element of the generic op- 
erator O is a linear combination of the matrix elements 
{^i \ O IV'i'). Such elements in turn are computed by ex- 



panding the CSFs on the SD basis [Eq. (7)]: 

The usage of the binary representation of SDs, as in 
Eq. ljA17p . allows for extremely efficient bit-per-bit oper- 
ations when evaluating matrix elements. Consider, e.g., 
the matrix element 

(^/IcltciiCciCdTll'.), (A21) 

which might represent a term in the Coulomb interaction. 
The ket 

ci^iccicd^ \^^) (A22) 

will differ from 1$^) at most in the occupancy of one or- 
bital with spin up and one with spin down. Translated 
into the binary representation, this implies that in either 
the spin- up or -down integers identifying the ket HA22|I . 
a couple of bits at most will be swapped with respect to 
their values being changed from (01) to (10) or vice 
versa (or no change at all). The matrix element (|A21|) 
will differ from zero only if the ket (|A22() is equal to |$/) , 
but for the sign: this is efficiently checked by performing 
a bit-per-bit xor operation between the same-spin inte- 
gers representing and |$/), respectively. Only if the 
number of resulting bits set to true (after the xor opera- 
tion per each spin) is equal to or 2, the matrix element 
(IA21|) may be non null. Eventually, the sign of (jA21|) 
is determined by counting how many bits set to 1 occur 
between those swapped. With those and similar tricks 
all types of matrix elements are easily evaluated. 

The above procedure is especially advantageous for 
post-processing)Qi^^i"i^^i^^i^'^i^^iS'^i^^i^'^i^°^ since the eva- 
lutation of various matrix elements is particularly trans- 
parent and simple once physical operators are expressed 
in their second quantized form. 



APPENDIX B: COULOMB MATRIX ELEMENTS 
OF FOCK-DARWIN STATES FOR THE 2D 
HARMONIC TRAP 

We report, for the sake of completeness, the explicit 
expressions of the Coulomb matrix elements Vabcd re- 
ferred to Fock-Darwin orbitals (Sec. IIIA), which have 
been used in the calculations of few-electron states in a 
2D harmonic trap presented in the th is paper (Sec. V). 
For a full derivation see Refs. Il3llll3l 

Since each Fock-Darwin orbital [Eq. (12)] is identified 
by the radial and azimuthal quantum numbers n and to, 
respectively, the generic Coulomb matrix element is iden- 
tified by eight indices: V„imi,„2m2,n3m3,n4m4- Among 
such indices, only seven are independent, in virtue of to- 
tal orbital angular momentum conservation. The explicit 
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expression is: 



has been used. 



Vr 



ni mi ,7127712 ,^37713 ,7147714 ^^7711+7712 ,7713+7714 



n 



E 



(-1) 



(4)j=0 



ii!j2!j3!j4! 



^ {rii + |mi|) 



1=1 



m - ji 



(4)^ = 

11(7) r(A/2 + i) 

X r([G- A + l]/2). (Bl) 



Here, r(^) denotes the Gamma fmiction and we use the 
foUowing conventions: 

(i) The shorthand 



rii 712 "3 "4 



E -EEEE 

(4)j = Jl=0j2=0j3=0j4=0 



(ii) The 7's are defined as: 



71 

72 
73 
74 



Jl + J4 

h + js 
h + js 
ji + k 



{\mi\+m{) /2- 
{\m2\+m2) 
{\m2\-7n2) 12- 
(|mi| -mi)/2- 



- (|to4| - 7714) /2, 

- (Imal - TO3) /2, 

- (Imal +m3)/2, 

- (|to4| + nii) /2. 



(iii) G and A are defined as: 



G = 71 + 72 + 73 + 74, 

A = ^1 +£2 +4+^4. 
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